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Abstract 



We describe a series of simulations of high temperature QCD with two fla- 
vors of Wilson quarks aimed at clarifying the nature of the high temperature 
phase found in current simulations. Most of our work is with four time slices, 
although we include some runs with six and eight time slices for comparison. 
In addition to the usual thermodynamic observables we study the quark mass 
defined by the divergence of the axial current and the quark propagator in the 
Landau gauge. We find that the sharpness of the Nt = 4 thermal transition 
has a maximum around k = 0.19 and = 4.8. 

12.38.Gc, ll.15.Ha 
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I. INTRODUCTION 



Lattice simulations are an important source of information on the behavior of quantum 
chromodynamics at high temperature. Most work has been done with Kogut-Susskind 
quarks because of the exact remnant of chiral symmetry. Since the exact chiral symmetry 
of Kogut-Susskind quarks is a U{1) symmetry, there is some question about how weU the 
results reproduce the real world with its SU (2) chiral symmetry. In the continuum limit the 
complete chiral symmetry is restored. However, in the continuum limit the results should be 
independent of the regularization used for the quarks. To test this it is important to study 
high temperature QCD with the other common form of lattice quarks, the Wilson quarks. 

The first simulations of high temperature QCD with two fiavors of Wilson quarks revealed 
a potential problem — for the values of for which most low temperature simulations 
were done, 4.5 < G/g^ < 5.7, the high temperature transition occurs at a value of quark 
hopping parameter k, for which the pion mass measured at zero temperature is quite large 
[1,2]. In other words, it is difficult to find a set of parameters for which the temperature is 
the critical temperature and the quark mass is small. Further work confirmed that the pion 
mass is large at the deconfinement transition for this range of [3,4]. (A recent study 

has concluded that for four time slices the chiral limit is reached at a very small value of 

Screening masses for color singlet sources show an approach to parity doubling in the 
high temperature phase similar to what is seen with Kogut-Susskind quarks [2,3]. Also, 
measurements of the pion mass show a shallow minimum at the high temperature transition 
[6]. 

Previous simulations with Wilson fermions have located Kj, the value of the hopping 
parameter at whch the high temperature crossover or phase transition occurs, as a function 
of 6/g^ for = A and 6. The critical value of the hopping parameter, k^, for which 
the pion mass vanishes at zero temperature has been located with somewhat less precision 
[1,2,6,7,3,4]. Some measurements of hadron masses been carried out on zero temperature 
lattices for values of k, and 6/g^ close to the Hf curve, allowing one to set a scale for the 
temperature, and to estimate Kc in the vacinity of the thermal transition [3,4,8]. 

In more recent work at = 6 we have observed coexistence of the low and high tem- 
perature phases over long simulation times, and we have extended these observations in the 
present project. The change in the plaquette across the transition is much larger than for the 
high temperature transition with Kogut-Susskind quarks [4]. This unexplained behavior, as 
well as work by Hasenfratz and DeGrand on the effect of heavy quarks [9], has led us to 
extend our work. 

This paper reports on a series of simulations with Wilson quarks at high temperature, 
in which we have studied a number of indicators for the nature of the phases. Using 8"^ X 4 
lattices, we have extended earlier studies of the location of the thermal transition or crossover 
to K, = 0.20, 0.21 and 0.22. In the ranges of 6/g^ and k, that have been studied earlier, we 
have done extensive simulations on 8 X 8 X 20 X 4 lattices, with additional work on 12"^ X 6, 
12^ X 24 X 6 and 8^ x 20 x 8 lattices. For one value of 6/g^ we made a series of runs on 
6x6x20x4 lattices to make sure that the effects we see are not due to the spatial size of the 
lattice. In this report we will concentrate on the results with four time slices. Simulations 
with Nf = 6 and 8 are still underway and will be described later. 
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We see a number of inexplicable effects. At large /3 and small n the crossover from the 
confined phase to the high temperature phase is smooth. Beginning at n) = (5.f , .f6) and 
extending down to about n) = (4.5f , .20) the crossover becomes abrupt, though probably 
not first order. A rapid crossover is seen in the plaquette, real part of the Polyakov loop, 
?/'?/', the entropy, and the quark mass derived from the axial current. For fi < 4.5, n > 0.20 
the transition once again becomes very smooth. 

Section 2 discusses the quantities we measured, and section 3 summarizes the simulations 
and the results. Conclusions are in section 4. 



II. MEASURED QUANTITIES 

In our simulations we have measured the expectation values of the Polyakov loop, the 
space-space and space-time plaquettes, the chiral condensate t/jt/j^ the entropy, screening 
masses for meson sources, the quark mass defined by the divergence of the axial current, 
and quark propagators in Landau gauge. 

The expectation value of the Polyakov loop, (P), is simply interpreted as exp(— Fg/T), 
where Fq is the free energy of a static test quark. With dynamical quarks (P) is always 
nonzero, but it increases dramatically at the high temperature transition. We also measured 
the space-space and space-time plaquettes, (Qss) and {^st}- In our normalization these are 
equal to three on a completely ordered lattice. 

The energy, pressure, entropy and t/jt/j with Wilson quarks are obtained by differentiating 
the partition function with respect to the temporal size, the spatial size, and the quark mass, 
respectively. Details are given in appendix A. We study the entropy to lowest order in g 
and t/jt/j^ using the formulae: 



1 inNf / I \ 



and 




d 



an 

4\ I 4 I 4 

where Sg and sj are the gluon and fermion entropies, respectively. 

We measured screening masses for meson sources with quantum numbers of the tt, cr, 
and ai. These measurements are a standard hadron spectrum calculation, except that the 
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propagation is in the z direction. We used a wall source covering the entire 2; = slice of the 
lattice, with the gauge fixed to a spatial Coulomb gauge which maximizes the traces of the 
x, y and t direction links. After blocking five to ten measurements together to minimize the 
autocorrelations, we fit all the propagators to a single exponential using the full covariance 
matrix of the propagator elements. 

A quark mass can be defined from the divergence of the axial current [10,11]. The basic 
relation is a current algebra relation 

• (V'75V'(0)V'757AtV'(a;)) = 2m,(V'75V'(0)V'75V'(a;)) (4) 

If we sum over slices, and measure distance in the z direction, this becomes: 

d - - 

Zl(V'75V'(0)V'7573V'(a^)) 

x,y,t 

= 2m, ^(V'75V'(0)V'75V'(a;)) (5) 

x,y,t 

We define PS{z) as the pion correlator with a point sink: 

PSiz) = (W/(0)V'75V'(^)) (6) 
and A[z) as the axial current correlator: 

A{z) = {w{oyhs7M^)) (7) 

where 1^(0) is the wall source at 2; = 0. At long distances both PS{z) and A[z) will fall off 
as exp(— m^2;). Therefore we perform a simultaneous fit to the two propagators on a lattice 
periodic in the z direction using three parameters, C, and m,, 

PS{z) = C sinh(m^) [exp [—nij^z) + exp [mj^[Nz — z))] (8) 

A[z) = C 2 niq [exp [—nij^z) — exp [mj^[Nz — z))] (9) 

The factor of sinh(m^) in Eq. (8) comes from using the lattice difference f[z-\-l)—f[z — l) for 
the derivative in Eq. (5). Note that PS{z) is periodic in z while A[z) is antiperiodic. We use 
the pointlike axial current t/j[z)^f5^f^t/j[z) rather than a point split current. These are quark 
masses in lattice units; to convert to continuum quark masses requres a lattice-to-continuum 
renormalization. See Ref. [17] for a discussion of this point. 

The quark propagator in the Landau gauge was also measured. This propagator has 
been studied with Kogut-Susskind quarks in Ref. [12]. We chose a source constant in the y 
direction and a 6 function in z and t with only the real part of the first Dirac component 
non-zero (in the Weyl basis we use). Because of the (5-function all possible momenta in x, 
z and t directions were excited. To distinguish among the different momenta we performed 
a Eourier transform of the propagator in x and t directions (taking into account that it has 
to have odd frequencies in t direction). This gives the propagation of the quark in the z 
direction as function of k^^ and kt, i.e. the dispersion relation of the screening propagator. In 
order to keep the amount of generated data at a reasonable level, the propagator was saved 
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only lor on-axis momentum values ol and kf. This enabled us to measure the on- axis 
dispersion relation ol the quark screening mass, in particular the screening mass difference 
ol the quark and light doublers. 

The form to which the spatial propagator is fitted is usually motivated by the form ol the 
tree propagator. We suppose that at large distances, each separate momentum component 
ol the spatial propagator resembles the corresponding tree quark form, but with its own 
renormalized quark mass, or in this case ol Wilson fermions, with a renormalized k. 

In momentum space the tree Wilson propagator is 



G{k) 



[1 - 2k cos(p^)] - i2k If, sin(p^ 



(10) 



[1 - 2k cos(p,)]2 + E, sm\p^) ' 
With our choice ol the source the first Dirac component ol the propagator, Gi, is real: 

[1 - 2/€ cos(p^)] 



[1 - 2k Eu cos(p^)]2 + 4/€2 E„ sin^iPf, 



Then, for non-zero z values 
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Gipi,P2,z,po) = ^exp[z27rA;2;/iv^]G'(pi,p2,^,Po) 



k=i 



;i - 6/€ + 4/€A)2 - ak{b + 1; 

8/€(l - 6/€ + iKAf 



cosh[ma(2; — Lz/2)] 
sinh(ma)sinh(maiv^/2) 



where 



At=l,2,0 



:i2) 



:i3i 



B = J2 sin^b^*]' 

At=l,2,0 



(14) 



and 



ma 



2 sinh ^ 



4:kW + (1 - Sk + AkAY 
\\ 8k{1 -Qk + AkA) 



(15) 



At zero momentum (on a low temperature lattice where the lowest Matsubara frequency is 
close to zero) this relation turns into 



ma = In 



1 — 6k 



2k 



The mass vanishes when /€ — )• 1/8 as expected. Inverting this for k gives 

1 



2 exp[ma] -|- 6 ' 

which, for small masses reduces to the naive relation. 



(16) 



:i7) 
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1 / X 

K= 18 

2ma + 8 ^ ^ 

that one expects looking at the terms of the Lagrangian. For large lattices the lowest doubler 
mass becomes 



madovLhier = lim In 



1 - 2k 



2k 



ln[3] = 1.09861. (19) 



For Kogut-Susskind fermions [12] the free propagator turns out to be a sum of two terms, 
having parts with an alternating sign in 2;-direction. For Wilson fermions, with our choice 
of the source, the propagator is a single exponential, or hyperbolic cosine, on a finite lattice. 

Furthermore, the sign of Gi at A; = changes at Kc = 1/8. Therefore, measuring the 
sign of the propagator can be used as an indicator of whether k, is effectively greater or less 
than Kc- 

One can infer from Fqs. 13-15 that the only effect of finite spatial lattice size is the 
discretisation of the momenta. For a given momentum all lattice sizes give the same value of 
the screening mass. For a smaller lattice, the range of allowed momenta is more restricted, 
of course. 




FIG. 1. The spatial screening mass at different spatial momenta for free Wilson fermions as a 
function of k for an 8^ X 20 X 4 lattice. At Kc = 1/8 the higher masses are for higher momenta. 

To be specific let us look at what happens with our lattice size: 8^ X 20 X 4. This is 
shown in figure 1. At the lowest momentum screening mass is at its minimum. If one 
increases k, the screening masses start to converge to a single value close to one at k, = 0.152. 
At this point the dispersion relation is fiat. 

The sign of the propagator with this source depends on the momentum. Generally, the 
value of K, at which the sign changes increases with the momentum. For zero momentum 
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it occurs at Hc] for the smallest nonzero momentum in our lattice size it takes place at 
n = (2 + \/2 — a/3)/11 = 0.1374. The amplitude for the doubler does not change sign in this 
kappa range. 

For our lattice size, inserting the appropriate momenta to Eq. (15) one obtains the 
following quark screening masses at k.^ = 1/8): 



k = 


(0,0,0,7r/4) 


ma = 


0.6610 


k = 


(7r/4,0,0,7r/4) 


ma = 


0.8906 


k = 


(7r/2,0,0,7r/4) 


ma = 


1.1171 


k = 


(37r/4,0,0,7r/4) 


ma = 


1.2149 


k = 


(7r,0,0,7r/4) 


ma = 


1.2411 



(20) 



For purely temporal momenta the free field screening mass is 

A; = (0,0,0,±7r/4) , ma = 0.6610 , 
A; = (0,0,0, ±37r/4) , ma = 1.0711 ^ 

Hence, the temporal doubler's screening mass is smaller than that of the spatial doubler. 



III. SIMULATIONS AND RESULTS 







6/^2 


K 


traj. 


ignore 


dt 


accept 


4 


8 


4.9 


0.180 


650 


100 


0.02 


0.88 


4 


8 


4.9 


0.181 


320 


100 


0.02 


0.94 


4 


8 


4.9 


0.182 


810 


100 


0.02 


0.90 


4 


8 


4.9 


0.1825 


824(c) 


100 


0.02 


0.86 


4 


8 


4.9 


0.1825 


780(h) 


100 


0.02 


0.88 


4 


8 


4.9 


0.183 


624(h) 


100 


0.02 


0.90 


4 


8 


4.9 


0.183 


810(c) 


600 


0.02 


0.87 


4 


8 


4.9 


0.184 


610 


100 


0.02 


0.93 


4 


8 


5.0 


0.173 


540 


100 


0.025 


0.86 


4 


8 


5.0 


0.175 


500 


100 


0.02 


0.92 


4 


8 


5.0 


0.177 


400 


100 


0.02 


0.90 


4 


8 


5.0 


0.178 


474 


100 


0.02 


0.94 


4 


8 


5.0 


0.180 


630 


100 


0.02 


0.92 


4 


8 


5.0 


0.182 


256 


100 


0.02 


0.86 


6 


8 


5.0 


0.175 


240 


80 


0.02 


0.94 


6 


8 


5.0 


0.180 


360 


60 


0.0167 


0.87 


8 


8 


5.0 


0.175 


350 


100 


0.0167 


0.91 




TABLE I. 


Table of runs 


at fixed with 


varying k. 


"(h)" and "(c)" 


indicate hot 


and cold 



starts. 

Simulations were run on the Intel iPSC/860 and Paragon, and the nCUBF-2 at the San 
Diego Supercomputer Center, on the Thinking Machines Corporation CM5 at the National 
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6/^2 


K 


traj. 


ignore 


dt 


accept 


4 


8 


5.1 


0.165 


1120 


100 


0.025 


0.89 


4 


8 


5.1 


0.167 


2660 


100 


0.025 


0.89 


4 


8 


5.1 


0.169 


460 


100 


0.025 


0.87 


4 


8 


5.1 


0.170 


700 


100 


0.025 


0.90 


4 


8 


5.1 


0.171 


980 


100 


0.025 


0.89 


4 


8 


5.1 


0.172 


3380 


100 


0.025 


0.87 


4 


8 


5.1 


0.173 


500 


100 


0.025 


0.85 


4 


8 


5.1 


0.175 


460 


100 


0.025 


0.86 


4 


8 


5.1 


0.177 


1800 


100 


0.025 


0.88 


4 


8 


5.1 


0.179 


660 


100 


0.025 


0.88 


4 


6 


5.1 


0.169 


860 


200 


0.0333 


0.82 


4 


6 


5.1 


0.170 


1000 


100 


0.0333 


0.85 


4 


6 


5.1 


0.171 


1360 


100 


0.0333 


0.81 


4 


6 


5.1 


0.172 


1120 


100 


0.0333 


0.83 


4 


6 


5.1 


0.175 


720 


100 


0.0333 


0.81 


8 


8 


5.1 


0.167 


512 


100 


0.025 


0.82 


8 


8 


5.1 


0.173 


279 


100 


0.02 


0.87 


8 


8 


5.1 


0.177 


440 


100 


0.02 


0.66 


4 


8 


5.3 


0.155 


2400 


100 


0.0333 


0.78 


4 


8 


5.3 


0.157 


660 


100 


0.0333 


0.76 


4 


8 


5.3 


0.158 


1239 


100 


0.0333 


0.89 


4 


8 


5.3 


0.159 


660 


100 


0.0333 


0.89 


4 


8 


5.3 


0.160 


1777 


100 


0.025 


0.88 


4 


8 


5.3 


0.161 


480 


100 


0.025 


0.88 


4 


8 


5.3 


0.162 


480 


100 


0.025 


0.83 


4 


8 


5.3 


0.163 


680 


100 


0.025 


0.87 


4 


8 


5.3 


0.164 


460 


100 


0.025 


0.90 


4 


8 


5.3 


0.165 


720 


100 


0.025 


0.88 


4 


8 


5.3 


0.166 


660 


100 


0.025 


0.87 


4 


8 


5.3 


0.167 


912 


100 


0.025 


0.86 


4 


8 


5.3 


0.168 


540 


100 


0.025 


0.89 


4t 


8 


5.3 


0.168 


840 


100 


0.025 


0.80 


4 


8 


5.3 


0.169 


380 


100 


0.025 


0.86 


4 


8 


5.3 


0.170 


380 


100 


0.025 


0.79 


4 


8 


5.3 


0.172 


440 


100 


0.025 


0.87 


6 


12 


5.3 


0.155 


320 


100 


0.0177 


0.88 


6 


12 


5.3 


0.160 


552 


60 


0.0177 


0.91 


6 


12 


5.3 


0.165 


666 


216 


0.0177 


0.85 


6 


12 


5.3 


0.166 


1403 


400 


0.0177 


0.84 


6 


12 


5.3 


0.167 


760 


302 


0.0177 


0.84 


6 


12 


5.3 


0.168 


603 


200 


0.0177 


0.85 



TABLE II. Table of runs at fixed with varying k. "(h)" and "(c)" indicate hot and cold 

starts. Most of the iV^ = 4 runs were on 8^ X 20 X 4 lattices. The run indicated with a | at 

K = 0.168, = 5.3 was done on a 8^ X 40 X 4 lattice. 
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' x,y,z 


6/^2 


K 


traj. 


ignore 


dt 


accept 


4 


8 


4.75 


0.19 


578 


100 


0.014286 


0.912(13) 


4 


8 


4.755 


0.19 


1504 


750 


0.014286 


0.960(6) 


4 


8 


4.76(c) 


0.19 


837 


500 


0.014286 


0.948(11) 


4 


8 


4.76(h) 


0.19 


362 


100 


0.014286 


0.962(12) 


4 


8 


4.32 


0.20 


156 


50 


0.02 


0.83(4) 


4 


8 


4.36 


0.20 


172 


50 


0.02 


0.83(3) 


4 


8 


4.40 


0.20 


188 


50 


0.02 


0.79(3) 


4 


8 


4.44 


0.20 


152 


50 


0.02 


0.77(4) 


4 


8 


4.48 


0.20 


368 


50 


0.02 


0.69(3) 


4 


8 


4.50 


0.20 


244 


100 


0.014286 


0.84(3) 


4 


8 


4.52 


0.20 


841 


150 


0.014286 


0.773(16) 


4 


8 


4.54 


0.20 


566 


150 


0.014286 


0.72(2) 


4 


8 


4.56 


0.20 


1324 


200 


0.014286 


0.941(7) 


4 


8 


4.60 


0.20 


365 


50 


0.014286 


0.937(14) 


4 


8 


4.64 


0.20 


244 


50 


0.02 


0.959(14) 


4 


8 


4.10 


0.21 


74 


50 


0.01 


0.79(8) 


4 


8 


4.20 


0.21 


267 


50 


0.005 


0.90(2) 


4 


8 


4.26 


0.21 


586 


100 


0.005 


0.85(2) 


4 


8 


4.28 


0.21 


478 


100 


0.005 ^ 0.0025 


0.82(2) 


4 


8 


4.30 


0.21 


454 


100 


0.005 


0.904(16) 


4 


8 


4.32 


0.21 


227 


50 


0.005 


0.94(2) 


4 


8 


4.34 


0.21 


259 


50 


0.007143 


0.943(15) 


4 


8 


4.36 


0.21 


281 


100 


0.0025 ^ 0.007143 


0.960(14) 


4 


8 


4.40 


0.21 


197 


50 


0.0025 ^ 0.01 


0.95(2) 


4 


8 


4.44 


0.21 


249 


50 


0.007143 ^ 0.01 


0.977(9) 


4 


8 


4.50 


0.21 


120 


50 


0.05 ^ 0.01 


0.94(3) 


4 


8 


3.80 


0.22 


59 


15 


0.001 


0.93(4) 


4 


8 


3.90 


0.22 


56 


30 


0.002 ^ 0.0004 


0.69(9) 


4 


8 


3.96 


0.22 


39 


25 


0.002 ^ 0.0005 


0.36(13) 


4 


8 


4.00 


0.22 


119 


80 


0.004 ^ 0.002 


0.50(8) 


4 


8 


4.04 


0.22 


161 


50 


0.004 


0.71(4) 


4 


8 


4.06 


0.22 


119 


40 


0.003333 


0.78(5) 


4 


8 


4.10 


0.22 


234 


50 


0.005 


0.86(3) 


4 


8 


4.20 


0.22 


90 


50 


0.007143 


0.98(3) 


4 


8 


4.30 


0.22 


58 


50 


0.005 ^ 0.007143 


1.00(0) 


4 


8 


4.40 


0.22 


90 


50 


0.01 


0.95(3) 


4 


8 


4.50 


0.22 


122 


50 


0.01 ^ 0.02 


0.96(2) 




TABLE III. Table of runs 


at fixed 


K with varying 


"(h)" "(c)" indicate hot 


and cold 


starts. The accept 


ance rate gives the average over 


all runs in 


the sample kept for measurement, 


whether or not dt was changint 


2; during 


the runs. 
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Center for Supercomputing Applications, and on a cluster of RS6000 workstations at the 
University of Utah. We used the hybrid Monte Carlo algorithm with two flavors of dynamical 
quarks in all our simulations [13]. The parameters of our runs are listed in tables I, II and 
III. 

For the 8^ X 20 X 4 runs we used trajectories with a length of one unit of simulation 
time and made measurements after every second trajectory. The step size for these runs 
ranged in the normalization of Ref. [14] from 0.033 for the largest 6/(/^ and smallest n to 
0.02 at the other extreme. Acceptance rates for these runs range from 70% to 90%, with an 
average over all the runs of 87%. For computation of the fermion force in the updating and 
the propagators in the measurements we used the conjugate gradient algorithm with even- 
odd ILU preconditioning [15]. The conjugate gradient residual, defined as M'^Mx — b / \b\ 

where M is the preconditioned matrix, b is the source vector and x is the solution vector, 
was 10~^. Runs were made at = 5.3, 5.1, 5.0 and 4.9 with = 4. At = 5.3 we 
also made a series of runs with = 6. At = 5.1 we ran two points with = 8 and 

at = 5.0 two points with Nt = 6. We also ran a series of simulations at 6/g^ = 5.1 on 

6^ X 20 X 4 lattices to verify that the spatial size of the lattice was not seriously affecting 
our results. At 6/g^ = 5.3 a series of short runs on 6"^ X 4 lattices was made for very large k. 

For reference we show a phase diagram for the relevant range of k, and 6/g^ in Fig. 2. 
Previous work showed that as k, increased from 0.16 to 0.19 along the = A high tem- 
perature crossover line the pion mass decreased, suggesting a closer approach to the high 
temperature transition in the chiral limit [3]. More recent work by Iwasaki et al., beginning 
from the 6/(/^ = limit, suggested that a high temperature transition for zero quark mass 
might be found at /€ 0.225 [16]. We have done a series of runs on 8"^ X 4 lattices in which 
we varied 6/(/^ at /€ = 0.20, 0.21 and 0.22 to extend the previous work. As expected, the 
number of conjugate gradient iterations required in the updating increases as n increases in 
this range, and the size of the possible updating time step decreases. Thus these runs have 
limited statistics. In Fig. 3 we show the plaquette and Polyakov loop as a function of 6/(/^ 
for the various values of n. Notice that the transition appears to be sharpest at /€ 0.19, 
becoming smoother for larger and smaller n. Fven in those cases where the transition is 
very abrupt, we do not see the sorts of metastability and tunneling characteristic of strong 
first order transitions. We do find cases where equilibration takes a long time. The worst 
case was in the run at 6/(/^ = 4.9 and n = 0.1825. In this case we have plotted two points, 
from hot and cold starts. These points are marked by arrows in Fig. 3. However, these two 
runs eventually converged to similar values, lying in between the values in the early parts 
of the runs. The time history of the Polyakov loop in these two runs is shown in Fig. 4. 

We now examine the 8^ X 20 X 4 runs in more detail. Figure 5 shows the real part of 
the Polyakov loop as a function of n for the different values of 6/(/^. For 6/(/^ = 5.3 we 
also include values for A^j = 6 to show how the transition point moves as N-t increases. For 
all of these values of 6/(/^ we see the expected sharp increase in the Polyakov loop at a 
values of n less than k^, where is the value at which the squared pion mass vanishes on 
a zero temperature lattice. We estimate Hc at these values of 6/(/^ from published values 
of He in Refs. [7] and [4] and a recent measurement at 6/(/^ = 5.3 by the HFMCGC group: 
/€c(5.3) = 0.16794 [17]. From a quadratic fit to these values, shown by a line in Fig. 2, we 
find /€c(6/5^) = 0.1687(2) at 5.3, 0.1795(4) at 5.1, 0.1861(12) at 5.0 and 0.1941(40) at 4.9. 

Although not a physical quantity, the number of conjugate gradient iterations used in 
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FIG. 2. Phase diagram showing estimates for the high temperature transition and k^. Circles 
represent the high temperature transition or crossover for Nt = 4, squares the high temperature 
transition for Nt = 6 and diamonds the zero temperature Kc- Previous work included in this figure 
is from Refs. [2], [6], [7], [3] and [4]. We show error bars where they are known. For series of 
runs done at fixed k the error bars are vertical, while for series done at fixed the bars are 

horizontal. Points coming from this work are shown in heavier symbols. The solid lines are fits to 
Kt for Nt = 4 and to Kc used in interpolating and extrapolating. 
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4.0 4.5 5.0 5.5 

6/g^ 




0.0 ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ — ^ ' ^ 
4.0 4.5 5.0 5.5 

6/g^ 

FIG. 3. The plaquette and Polyakov loop as a function of 6/^^ for various values of k. The 
diamonds are previous results of Ref. [3] for k = 0.12, 0.14, 0.16, 0.17, 0.18 and 0.19. For k = 0.12 
and 0.14 data from long runs as well as some data from short runs collected while generating 
hysteresis loops is shown. The octagons at k = 0.20, 0.21 and 0.22 are new results from 8"^ X 4 
lattices. The squares come from runs on 8^ X 20 X 4 lattices. These runs were done at fixed values 
of with varying k. They have been mapped onto this figure by fitting the 6/^|, Kt line (with 
a fit shown as a line in Fig. 2), and moving the points in the k, plane parallel to this line. 

Specifically, we plot the points at ^jglffective = ^/Srun ~ °'^g{'^*^ (Krun - Kt). The fit for at 
6/^2 ^ 5 5 and 4.9 is 0.1579, 0.1713, 0.1772 and 0.1827 respectively. 
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FIG. 4. Time history of the real part of the Polyakov loop for runs with hot and cold starts at 
6/^2 ^ 4 9 ^ ^ 0.1825. 

solving M^Mx = b indicates how singular M is on the average. This quantity has been used 
as a probe of the physics in Ref. [16] In Fig. 6 we show the average number of conjugate 
gradient iterations used in an updating step, where a linear extrapolation of the last two 
time steps was used to produce a starting guess for the solution vector. For = 5.3 and 
Nf = 4: there is very little effect on the number of iterations at k.^. As is decreased for 
Nf = 4: there is an increasingly sharp peak in the number of iterations at k.^. Notice also the 
sharp peak in the = 6 results for = 5.3. 

Figure 7 shows the average plaquette in these runs. Our normalization is such that the 
plaquette is three for a lattice of unit matrices. The plaquette also shows a sharp rise as 
the high temperature crossover is passed. Notice that for = 5.3 we have results for 

Nf = 4: and 6 showing that this increase is in fact due to the time size of the lattice, or the 
temperature. 

The chiral condensate t/jt/j is less useful for Wilson quarks than for Kogut-Susskind quarks, 
since it does not go to zero in the high temperature phase without difficult subtractions. 
Nevertheless we plot it in Fig. 8. There is a clear drop in t/jt/j as the high temperature 
transition is crossed. This drop increases dramatically as 6/g^ decreases. 

Perhaps the most physically relevant observable is the entropy. In Fig. 9 we plot T X s 
in units of a""*. To give an idea of the normalization of this graph, for eight gluons in 
free field theory on an 8 X 8 X 20 X 4 lattice the gluon entropy would be Tsgiuejree = 
0.0400""*, while for two fiavors of free Wilson quarks at k, = Kc = 0.125 the entropy would be 
Tsquarkjree = 0.1250""*. The cffects of the lattice spacing and spatial size are very large here; 
in the continuum with infinite spatial extent these numbers are 0.027 and 0.036 respectively. 
Strangely, when we divide the entropy into gauge and fermion parts as in Eqs. 2 and 3 we 
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FIG. 5. Expectation value of the Polyakov loop as a function of k for the various values of 
Results are shown for Nt = 4 for = 5.3, 5.1, 5.0 and 4.9 (octagons). For = 5.3 

we also show results for Nt = 6 (diamonds). The crosses along the = 5.1 line are results 

on a 6^ X 24 X 4 lattice at 6/g^ = 5.1, to show that the spatial size of the lattice is not greatly 
affecting the results. The dotted symbols extending the 6/g^ = 5.3 line are short runs on a 6"^ X 4 
lattice, showing that the behavior is smooth out to very large k. The vertical lines mark the zero 
temperature Kc for 6/g^ = 5.2, 5.1 and 5.0 respectively. (For 6/g^ = 4.9, Kc(T = 0) ~ 0.194.) 
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FIG. 6. Conjugate gradient iterations for updating step, as a function of k for the various values 
of 6/^^. Again the diamonds are Nt = 6 results at = 5.3. 
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FIG. 7. Expectation value of the plaquette as a function of k for the various values of 
Here we included values for larger Nt to emphasize the effect of the temperature. The dotted 
symbols for = 5.3 are short runs on a 6"^ X 4 lattice extending to k far beyond the zero 

temperature Kc shown by the vertical line. 
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Nt=4,6,B (0,0, x) 




0.15 0.16 0.17 0.18 0.19 

K 

FIG. 8. Expectation value of 'ipip as a function of k for the various values of 6/^^. The dotted 
symbols for 6/^^ = 5.3 are short runs on a 6"^ X 4 lattice extending to k far beyond the zero 
temperature Kc shown by the vertical line. 
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FIG. 9. Entropy (actually T X sa^) as a function of k for the various values of 6/^^. Again we 
show the 6^ X 24 X 4 results for comparison. 
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find that the gauge entropy is comparable to the fermion entropy instead of much smaller 
as would be the case with free fields on a lattice of this size. 
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FIG. 10. Quark mass from the axial current as a function of k for the various values of 
Points marked with question marks indicate runs where we were unable to get consistent fits as a 
function of distance. The plus signs on the niq = line are the zero temperature Kc for = 5.3, 
5.1 and 5.0. 

The quark mass defined by the divergence of the axial pion propagator is plotted in 
Fig. 10. When this quark mass was small we had great difficulty in getting good fits to 
the forms in Eqs. 8 and 9. This is expected, because when the quark mass is small the 
amplitude for the propagator A[z) is very small. Additionally, there is a tendency for the 
effective quark mass, or the quark mass coming from a fit over a short distance range, to 
increase with distance from the source. In cases where we were unable to get a fit with a 
satisfactory or where the quark mass was not convincingly independent of distance, we 
plot the point with a question mark in Fig. 10. To pursue this further we ran one of the 
difficult points, = 5.3 and k, = 0.168, on a 8^ x 40 x 4 lattice, allowing us to measure 

the ratio out to a distance of twenty. Figure 11 summarizes the results. In this figure we 
show the effective pion mass obtained from PS{z) and A[z) by fitting two two successive 
distances, and the quark mass obtained from simultaneously fitting both propagators at the 
two successive distances (a one degree of freedom fit). Unfortunately, in all other cases the 
lattice was only twenty sites long and we have to draw conclusions from distances less than 
ten. In Fig. 10 we see that when the Nf = 4: lattice enters the high temperature regime 
the pointlike axial current quark mass no longer agrees with the low temperature lattices 
[Nf = 6 and 8). The plusses at niq = in Fig. 10 are estimates for the zero temperature 
Kc- The axial current quark masses go through zero at k, less than the zero temperature 
Kc- When the axial current quark mass vanishes, the system is in the high temperature 
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FIG. 11. Pion effective screening masses from PS{z) (circles) and from A{z) (squares), and the 
effective quark screening mass from their ratio. The results are from an 8^ X 40 X 4 lattice with 
Qjc/ = 5.3 and k = 0.168. 



phase for /3 > 5.0, while at /3 = 4.9 kt apears to coincide with the point where the axial 
current quark mass vanishes, within experimental uncertainty. Note however, that the pion 
screening mass in the confinement phase is still nonzero at the transition point at /3 = 4.9. 

In Fig. 12 we show the squared pion screening masses in these runs. Again we see an 
increasingly sharp dip at Hf as decreases and k, increases. The appearance of the cusp at 
/3 = 5.1 coincides with the beginning of the region where the transition is abrupt. Screening 
masses for the tt, /?, a and ai mesons are shown in Figs. 13, 14 and 15. In all cases we see the 
screening masses coming together as the high temperature transition is crossed. However, 
we do not see any indication that the tt — a oi p — ai splittings in the high temperature 
regime are decreasing as Q/g^ decreases. Although the smaller pion masses in the cold regime 
suggest that chiral symmetry is being approached as we move toward smaller along 
the Kf line, we do not see this trend in the high temperature screening masses. Also notice 
that there are nonzero splittings between the parity partners at the points where the axial 
current quark mass is zero. Thus the vanishing of this quark mass is not an indicator for 
complete chiral symmetry restoration in the system. 

To investigate the contributions of the doublers to thermodynamic quantities such as 
the entropy we measured the effective masses from the quark propagator in Landau gauge 
at a few values of k, and G/g^. We find that fitting the quark screening propagators is more 
difficult than fitting the meson propagators. In part this is because the quark propagators 
fiuctuate more from configuration to configuration. There also seems to be a systematic 
trend toward larger effective quark masses at larger distances. With these caveats, the 
masses of the physical quark and the lightest doublers are given in Table IV. The fits were 
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FIG. 12. Pion screening mass squared as a function of k for the various values of 6/^^. Again, 
the circles are for Nf = 4, the diamonds for iV^ = 6 and the crosses for Nf = 8. The bursts are zero 
temperature pion masses from the HEMCGC collaboration at 6/^^ = 5.3. 



0=7T, □ = CT, 0=p, -''•i = S.^ 




0.155 0.160 0.165 0.170 



FIG. 13. Meson screening masses for 6/^^ = 5.30. The points connected by solid lines are for 
iVj = 4 and the points connected by dashed lines for Nf = 6. 
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FIG. 15. Meson screening masses for Qjg^ = 4.90 at Nf = 4. The two points at k = 0.1825 are 
from cold and hot starts. 
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selected by choosing the largest fit range that gives an acceptable confidence level. The 
ranges and confidence levels are also given in Table IV. 
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TABLE IV. The screening masses 


for the quark and the lightest 


doublers. 




is the 



difference of the spatial (temporal) doubter screening mass to the quark screening mass at the 
lowest momenta. The sign is for G(kt = vr/4). The fits were done simultaneously to all three 
propagators taking into account the cross correlations. The confidence level q and the range of 
each fit is also displayed. 



IV. CONCLUSIONS 

The most naive expectation regarding the thermodynamics of two fiavors of Wilson 
quarks at fixed Nt is that there would be a line in the /€, /3 plane at which a confinement- 
deconfinement transition occurs, that the transition would be smooth (crossover or second 
order), that the pion mass would smoothly decrease along that line, and that at some point, 
possibly corresponding to the point where the transition line crossed the zero temperature 
Kc — /3 line, the pion mass would go to zero. At that point one would have a finite temper- 
ature confinement-deconfinement or chirally restoring transition analogous to that seen in 
staggered fermions. Simple arguments [9] would put this point around /3 = 5.0 at = 4. 

These naive expectations are not borne out by the data. The chiral limit is reached 
at a very small /3 value if it is reached at all. However, near /3 = 5.0 A^t = 4 Wilson 
thermodynamics displays a number of features which have no analogs in staggered fermion 
systems. The transition becomes very sharp, though not first order as far as we can tell. A 
cusp in the pion screening mass appears as one crosses from the confined to the deconfined 
phase. The axial vector quark mass becomes strongly dependent at this point and for 
small Nf does not go to zero at its zero temperature k, value (at fixed /3). The sharp transition 
persists down to /3 = 4.5, k = 0.20 or so, at which point it is once again becomes smooth. 
As far as we can tell, the zero temperature n = Hc point plays no role in any A^t = 4 effects 
we have observed. 

It is tempting to speculate that the crossover line in the /€, 6/(/^ plane is close to some 
phase boundary where the transition is steepest. We are currently exploring this region with 
Nt = 6, where preliminary results indicate a change in the nature of the high temperature 
transition around this value of n. 

Indicators for the nature of the high temperature phase give a somewhat mixed picture. 
It is clear from the meson screening masses and from that chiral symmetry is at least 
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partially restored at high temperature. While the axial current quark mass goes to zero 
the TT — cr and p — ai splittings in the screening masses remain nonzero. Quark propagators 
in the Landau gauge suggest a large constituent quark mass at the transition, at least for 
= 5.3 and 5.1. This is consistent with earlier work [3] where at the = A crossover 
point near these k.) values the zero temperature pion was found to be quite heavy. 

Notice that the series of runs at = 5.3 extends to k, significantly larger than the zero 
temperature and there is no noticeable effect on any of the measured quantities when 
this Kc is crossed. (In fact, we have done short runs on 8"^ X 4 lattices for k, as large as 0.19 
at = 5.3 and seen no effects.) The sign of the propagator of the zero momentum quark, 
shown in table IV, is consistent with the sign of the axial current quark mass. Both of these 
quantities are behaving in the way one would expect in a free field theory at k, > n^- 
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Expressions for the energy, pressure and are found by differentiating the partition 
function with respect to 1 /T, volume and quark mass, respectively. Eirst, we write the action 
with adjustable lattice spacings in all directions. Introducing dimensionless parameters a^, 
we write the lattice spacing in the p direction as = aa^. Clearly this is redundant, since 
we have five parameters, a and to specify four lattice spacings, but it is convenient and 
symmetric. In the conventional notation of Karsch, = ai/ctf, where all the spatial a's are 
the same. When we are done taking derivatives, all the will be set to one. 

The partition function is 
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where the gauge action is 
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^. = EE4^^^^^T^°.^ (23) 

where D^j, is the plaquette in the iJ.v plane normahzed to three for unit matrices. We aUow 
a different gauge couphng g^i, in each plane. The fermion action is 

Sf = ^TTlogM^M (24) 

where 

M = 1-kJ2—P. (25) 

where 

P, = {1 + 7,W,{x)Sy,,+^ + (f - j.Wlix - j,)Sy,,_-, (26) 

The in the coefficient ofp^ takes care of the dimensional scaling of the first derivative. 
Notice that we have made a somewhat arbitrary choice in M when we scaled the irrelevant 
second derivative part with in the same way that we scaled the first derivative part. The 
in the coefficient must be adjusted to get correlation functions to be Euclidean invariant. 
Its role is similar to the Karsch coefficents C^j and Cr in the gauge action. Presumably has 
a power series expansion in g just as C^j and Cr- Once again we have more parameters than 
we need: four and k, for four directions. This parameterization is convenient because it 
includes the customary k, and, later, k.^ as parameters. We can fix the ambiguity well enough 
for our purposes by requiring that = I when all the are equal. In other words, if all 
directions are scaled by the same factor the only thing that changes is k. 

Let Kc be the value of k, at which the pion mass and quark mass vanish, at least on an 
infinite lattice. Following free field theory, we introduce a quark mass 

2ma = - (27) 

so that 

M= . \ {2ma + -T,^]D„ \ (28) 



^ + 2ma \ V "a* 



Here k,~^ will be a function of the couplings (/^j, and the scale factors a^. In free field theory, 
K,~^ = 2 Y^fM ct^^- (e^t = 1 in free field theory.) 

We find the energy, pressure and t/jt/j by differentiating the partition function: 

(29) 

\V,m constant V / 



V d/3 



P ^ ()\^ 1/^'"^ constant V / 



V^V^ — 7, \l3,V constant [■Ji-) 

pV dm 
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Here V is the volume, V = Hi Niai = N^a^. fi is the inverse temperature, fi = 
aNtdf Here the energy and pressure derivatives are taken with m constant, rather than 
with n constant. This is because Hc depends on the a^, so that if we distort the lattice while 
holding n fixed, the quark mass, and every physical mass, will vary sharply. 

We change the temperature and volume using 

- = — — (32) 

and 

' ' ' (33) 



dV m^a^al da. 



Alternatively, it may be easier to vary the volume by varying only one of the spatial lattice 
spacings 

— = 1 d 

dV N^a^a^ay da^ 

The gauge energy and pressure are standard [18-20]. Following Ref. [18], we define two 
derivatives. 

= (35) 

where A is not one of fi oi p and 

C\ = (36) 

where A is one of fi oi p. 

Because stretching both the time and space directions is equivalent to changing the 
lattice spacing, C^j and Cr are related to the beta function. 

^] 9-J = -2C\ - 2C\ = (37) 
\^ dlog(a) 

The contributions to the energy and pressure from Sg and Sj add. Doing the differenti- 
ation, and then setting the to one, the gluon energy is 

= ^ E i°st - o,,) + ean,, + 6C.n,,^ (38) 

Here and are the space-space and space-time plaquettes, again normalized to three 
for a lattice of unit matrices. 

For the gluon pressure, we find 

^5«'= 1^ E ( 4 (°«* - °-) - i^ss + 2n.t) - 2Cr (2n,, + a,,)) (39) 

^^.S^^t X \9 / 
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We also consider the linear combination e + the entropy. 



4 -'- 4 I 4 

aNt 

= l^i:(g+4,CV-a)) ,□.,-□„,) ,40, 

The entropy is obviously zero at T = 0. 

Just as the gauge couplings vary with the lattice spacings, k,~^ and vary with the 
lattice spacings as we try to hold m fixed. There is an explicit dependence of k,~^ on m plus 
a dependence of k,~^ on </, where g is varying with the a^. Now it clearly doesn't matter 
which direction we stretch the lattice, since k,~^ is defined on the infinite lattice, so 

oat oa^ 

Therefore 

There are two independent derivatives of the e^, analogous to and Cr- Define 

Br = ^ (43) 

and 

B. = ^. l^ + v (44) 

To compute we can just set and to one at the beginning. 

1 Nf I d 



a^tPtP= — — ^ ( Trlog M^M 



mNf 2 \d 



ma 



NfNt 2 \ MW \dma dma J / 

= ^^/Trf-L^ + ±^)\ ,45) 

NfNt 2 \ dma M dmaj I ^ ' 

The two parts are complex conjugates, so keep only one and take twice the real part. Using 
Eq. (28) 



1 iV^_ 1 / -2 . 1 



-2Re Tr 2ma + K:'-j:Pl]+2 



NfNt 2 \ ikft V(/€-i + 2ma)2 V = ^ " J \K-^ + 2ma 



1 4:KN 



f 



Re/Trf-^-1)) (46) 



mNf 2 \ VMt 
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where in the last step we used + 2ma = k. Looking at the derivation shows that 

the 1 in Eq. (46) comes from differentiating the + 2ma outside the parentheses in 

Eq. (28). Had we taken the fermion matrix to be 



M 



(47) 



this term would be absent. Since this latter form is closer to the usual continuum Lagrangian, 
we prefer 



NfNt 2 \ MU 



(48) 



as our expression for t/jt/j. 

Now for the fermion energy, differentiate Eq. (25) and then set and to one: 



1 Nf / 1 dM^\ 

^ /2ReTr-^- — ) 



NfNt 2 
-2 Nj 



Re/^ 
N^Nt 2 XMt 

f 



-f 



' K ^ + 2ma] 



2ma + - 



dat 



KNf 



NMt 2 



Re 



'Ok: 



dat 



(49) 



(50) 



where in the last step we used 1/ k,^ ^ + 2ma = n. 
Eor the fermion pressure: 



^ /2ReTr^-^ — ) 



3iV3iVt 2 
2 Nf 



Mt da, I 



Re/^ 

3iV3iVt 2 XMt 

f 



-f 



[k^ ^ + 2ma] 



2ma + - 



Ok: 



da. 



^ + 2ma \ ~~ oa, ^ oa. 



-2 kN 



(51) 



smN 2 



Re ( 1^ (l - ^) + ^ {3B.Pl + i2B. + 5.) - J:PI 



(52) 



Just as for t/jt/j^ the f term in ej and comes from differentiating the overall factor of 
K,~^ . It can be included or not, as desired. It will cancel when the finite parts of the energy 
and pressure are calculated by subtracting the zero temperature result from the nonzero 
temperature result. However, the j^-^^ term will not cancel out, since is temperature 
dependent. In practice (Tr I/Aft) is fairly close to (Tr I) = 4 X 3, so numerically it may be 
best to leave the I in. Then we would use 
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to express the energy and pressure just in terms of the expectation values of the spatial and 
temporal components oip'^ . 

Much of the difficulty cancels out if we look at the entropy, or sum of energy and pressure: 

f 

- = eja' + PfO* 

t 

'1 + {B.-B^))(iP,-\Y.P^) (54) 
The relation 

Oag oat 

resulted in all the derivative terms cancelling. (Remember, by d j da^ we mean d/da^ + 
d I day + d I daz — vary all the spatial ai together.) As usual, the entropy is obviously zero 
at r = 0, where |), =p Q. Since no zero temperature subtraction is required for the entropy, 
the terms involving i?^ and B„ will be higher order in than the "f" term, and we have 
neglected them in Eq. (3). 

g -1 

Obviously, the big problem in getting the energy and pressure separately is to find 

g -1 . . . ''i 

and and B„. Erom Euclidean invariance, is independent of ji. The variable k~ 

depends on the a^ in two ways. Eirst, there is an "explicit" dependence. Erom examining 
the fermion matrix, Eq. (25), we see that if all the a^ are scaled together with g held fixed. 
He is proportional to a. Thus 

explzczt 

Secondly, there is an "implicit" dependence of k,~^ on coming from the fact that Kc 
depends on g~^^ and we adjust the as we adjust the a^. Again, Euclidean invariance 
says 

dK^ = IJ^^ ^ (57) 

da. ■ 4 91og(a) 

SO only the beta function appears. 



dan ■ 4 91og(a) 

f dK.~^ d6/g'^ 



(58) 



idQ/g'^ dlog{a) 

We could proceed by estimating dn~^ jd&jg^ from our (i&it&i Hit various values of 6/(/^, or 
from correlations of the hadron propagators with the plaquette. Similarly, we could take a 
beta function either from perturbation theory or from some set of lattice simulations. We 
will not solve this problem here, so we will only quote the entropy rather than the energy 
and pressure separately. 



28 



REFERENCES 

[1] M. Fukugita, S. Ohta and A. Ukawa, Phys. Rev. Lett. 57, 1974, 1986. 
[2] A. Ukawa, NucL Phys. B (Proc. SuppL), 9, 463, 1989. 

[3] The HEMCGC Collaboration, K.M. Bitar, et a/., Phys. Rev. D, 43, 2396, 1991. 

[4] C. Bernard, T.A. DeGrand, C.DeTar, S. Gottlieb, A. Krasnitz, M.C. Ogilvie, R.L. Sugar 

and D. Toussaint, Phys. Rev. D, 46, 4741, 1992. 
[5] Y. Iwasaki, L. Kanaya, S. Sakai and T. Yoshie, to be published in Nucl. Phys. B [Proc. 

SuppL] . 

[6] K. Bitar, A.D. Kennedy, and P. Rossi, Phys. Lett., B234, 333, 1990. 
[7] R. Gupta, A. Patel, C. Baillie, G. Guralnik, G. Kilcup, and S. Sharpe, Phys. Rev. D 
40, 2072, 1989. 

[8] K. Bitar et a/., Nucl. Phys. B (Proc. suppl.) 30, 401, 1993. 
[9] A. Hasenfratz and T.A. DeGrand, Colorado preprint HEP-311, 1993. 
[10] M. Bochicchio, L. Maiani, G. MartineUi, G. Rossi and M. Testa, Nucl. Phys. B262, 331, 

1985; L. Maiani and G. MartineUi, Phys. Lett. 178B, 265, 1986. 
[11] S. Itoh, Y. Iwasaki, Y. Oyanagi and T. Yoshie, Nucl. Phys. B274, 33, 1986; Y. Iwasaki, 

L. Kanaya, S. Sakai and T. Yoshie, Phys. Rev. Lett. 67, 1494, 1991. 
[12] G. Boyd, Sourendu Gupta and F. Karsch, Nucl. Phys. B385 (1992) 481. 
[13] S. Duane, A.D. Kennedy, B.J. Pendleton, and D. Roweth, Phys. Lett., 195B, 216, 1987. 
[14] . Gottlieb, W. Liu, D. Toussaint, R.L. Renken and R.L. Sugar, Phys. Rev. D, 35, 2531, 
1987. 

[15] T. DeGrand, Comp. Phys. Comm., 52, 161, 1988; R. Gupta, A. Patel, C. Baillie, G. 

Guralnik, G. Kilcup, and S. Sharpe, Phys. Rev., D, 40, 2072, 1989; T. DeGrand and 

P. Rossi, Comp. Phys. Comm., 60, 211, 1990. 
[16] Y. Iwasaki, K. Kanaya, S. Sakai and T. Yoshie, Nucl. Phys. B [Proc. SuppL] 26, 311, 

1992, Phys. Rev. Lett. 69, 21, 1992. 
[17] The HEMCGC Collaboration, K.M. Bitar, et al., in preparation. 
[18] E. Karsch, Nucl. Phys. B205, 285, 1982. 

[19] J. Engels, E. Karsch, I. Montvay and H. Satz, Phys. Lett. lOlB, 89, 1981, Nucl. Phys. 

B205, 545, 1982; T. Celik, J. Engels and H. Satz, Phys. Lett. 129B, 323, 1983. 
[20] R. Trinchero, Nucl. Phys. B227, 61, 1983. 



29 



